*
* Copyright (c) 2020 NVI, Inc.
*
* This file is part of VLBI Field System
* (see http://github.com/nvi-inc/fs).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
      SUBROUTINE XAT (I0,E0,W0,O0,M0,A0,N0,I0Y,T0D,
     .                IMJD,UT,ALAT,ALON,HA,DEC,RANGE)
      IMPLICIT real*8 (A-H,J-Z)
      real*8 I0
      INTEGER K,JULDA
      include "../skdrincl/constants.ftni"
C
C     I0     :INCLINATION AT TIME T0 (DEG)
C     E0     :ECCENTRICITY (UNITLESS)
C     W0     :ARGUMENT OF PERIGEE (DEG)
C     O0     :R.A. OF ASC NODE (DEG)
C     M0     :MEAN ANOMALY (DEG)
C     A0     :SEMIMAJOR AXIS (KM)
C     N0     :MEAN MOTION (REV/SIDEREAL DAY)
C     I0Y    :EPOCH YEAR
C     T0D    :EPOCH DAY AND FRACTION
C
C     A0 AND N0 ARE MUTUALLY EXCLUSIVE. USE ONE AND SET THE OTHER TO
C     ZERO (BUT USE A VARIABLE, NOT LITERAL "0.D0", SINCE THE
C     SUBROUTINE RETURNS THE VALUE)
C
C     IMJD   :MODIFIED JULIAN DAY FOR CALCULATION.
C     UT     :UNIVERSAL TIME ON IMJD FOR CALCULATION.
C     ALAT   :LATITUDE OF OBSERVER (RADIANS)
C     ALON   :LONGITUDE OF OBSERVER (RADIANS)
C     HA     :APPARENT LOCAL HOUR ANGLE
C     DEC    :APPARENT DECLINATION
C     RANGE  :DISTANCE (KM)
C
C     TRANSLATION OF TAC AMSAT PROGRAM  BY CAK 830408
C     MODIFIED FOR USE BY SKED          BY WEH 830523
!     JMGipson 2004May14 constants moved outside of this routine.
C
C
C     GET TIME SINCE EPOCH
C
      I0DY=T0D
      I0YR=I0Y-1900
      IMJD0=JULDA(1,I0DY,I0YR)
C
      TDIF=IMJD-IMJD0
      TDIF=TDIF+(UT/secperday)
      TDIF=TDIF-(T0D-DBLE(FLOAT(I0DY)))
C
C   FIND SIDERAL TIME FOR EPOCH OF CALCULATION
C
      CALL SIDTM(IMJD,GSTM,FRACT)
      GSTM=(GSTM/(2.D0*PI))+SRATIO*(UT/secperday)
C
C
C     EVALUATE OBSERVERS GEOCENTRIC COORDINATES
C
C
C     X-TOWARD GREENWICH
C     Y-IN R.H.SENSE
C     Z-TOWARD N. POLE
C
      LAT=ALAT
      LON=ALON
      S9=DSIN(LAT)
      C9=DCOS(LAT)
      C8=DCOS(LON)
      S8=DSIN(-LON)
      R9=ERAD*(1.D0-(EFLAT/2.D0)+(EFLAT/2.D0)*DCOS(2.D0*LAT))
      L8=DATAN((1.-EFLAT)**2*S9/C9)
      X=R9*DCOS(L8)*C8
      Y=R9*DCOS(L8)*S8
      Z=R9*DSIN(L8)
C
C
C     UPDATE ELEMENTS TO CURRENT EPOCH.
C
C     FIRST CALCULATE MEAN MOTION FROM SEMIMAJOR AXIS OR VICE VERSA,
C     DEPENDING ON WHAT WAS SUPPLIED
C
      IF (N0.GT.0.1D0) A0=((GME/(N0**2))**0.333333333333333D0)
      IF (N0.LT.0.1D0) N0=DSQRT(GME/(A0**3))
      E2=1.D0-E0*E0
      E1=DSQRT(E2)
C
C     Q0    :INITIAL ORBIT PHASE
C
      Q0=M0/360.D0
C
C     MAKE ADJUSTMENT TO THE NODE TO ACCOUNT FOR PRECESSION DUE TO THE
C     EARTHS OBLATENESS (APPROXIMATELY)
C
      K2=9.95D0*((ERAD/A0)**3.5)/(E2**2)
C
C     PRECESSION OF THE NODE CALCULATED ELSEWHERE
C
      K2=0.
C
C     SIN AND COS OF ANGULAR ELEMENTS OF DATE
C
      FACT=PI/180.D0
      ARG=I0*FACT
      S1=DSIN(ARG)
      C1=DCOS(ARG)
      O=O0-TDIF*K2*C1
      ARG=O*FACT
      S0=DSIN(ARG)
      C0=DCOS(ARG)
      W=W0+TDIF*K2*(2.5D0*(C1**2)-0.5D0)
      ARG=W*FACT
      S2=DSIN(ARG)
      C2=DCOS(ARG)
C
C     ELEMENTS OF ROTATION MATRIX
C
      C11= C2*C0  -  S2*S0*C1
      C12=-S2*C0  -  C2*S0*C1
      C21= C2*S0  +  S2*C0*C1
      C22=-S2*S0  +  C2*C0*C1
      C31= S2*S1
      C32= C2*S1
C
C     MEAN ANOMALY OF DATE
C
C     K= PERIGEE PASSAGE COUNTER
C
      Q=N0*TDIF+Q0
      K=Q
      M=(Q-K)*2.D0*PI
C
C     ITERATE TO SOLVE KEPLER'S EQUATION FOR THE ECCENTRIC ANOMALY
C     GIVEN THE MEAN ANOMALY AND CALCULATE GEOCENTRIC COORDINATES
C     OF THE SATELLITE
C
      E=M+E0*DSIN(M)+0.5D0*(E0*E0)*DSIN(2.D0*M)
      ICT=0
2     CONTINUE
      ICT=ICT+1
3     CONTINUE
      S3=DSIN(E)
      C3=DCOS(E)
      R3=1.D0-E0*C3
      M1=E-E0*S3
      M5=M1-M
      IF (DABS(M5).LT.1.D-6) GO TO 1
      E=E-M5/R3
      GO TO 2
1     CONTINUE
C
C     GET SATELLITE POSITION IN ORBITAL PLANE
C
      X0=A0*(C3-E0)
      Y0=A0*E1*S3
      R =A0*R3
C
      LEN0=DSQRT(X0*X0+Y0*Y0)
C
C     ROTATE TO CELESTIAL COORDINATES
C
      X1=X0*C11+Y0*C12
      Y1=X0*C21+Y0*C22
      Z1=X0*C31+Y0*C32
C
C     CHECK ROTATION
C
C
C     ROTATE SO AXES COINCIDE WITH GEOCENTRIC AXES
C
      G7=GSTM
      IG7=G7
      G7=(G7-IG7)*2.D0*PI
      S7=-DSIN(G7)
      C7=DCOS(G7)
      XS=X1*C7-Y1*S7
      YS=X1*S7+Y1*C7
      ZS=Z1
C
C
      LEN2=DSQRT(XS*XS+YS*YS+ZS*ZS)
      LEN1=DSQRT(X1*X1+Y1*Y1+Z1*Z1)
      LEN3=DSQRT(C7*C7+S7*S7)
C
C     GET VECTOR FROM OBSERVER TO SATELLITE
C
      X5=XS-X
      Y5=YS-Y
      Z5=ZS-Z
      R5=DSQRT(X5*X5+Y5*Y5+Z5*Z5)
      RANGE=R5
C
C
C   GET RA AND DECK FROM GEOCENTRIC COORDINATES
C
      C5=Z5/R5
      S5=DSQRT(1.D0-C5*C5)
      DEC=.5D0*PI-DATAN2(S5,C5)
      HA=-DATAN2(Y5,X5)-LON
C
      RETURN
      END
